# This code replicates Table 7
rm(list=ls())
setwd("")
library(foreign)
data <- read.dta("data.dta")

data$religion_adventist <- ifelse(data$church_membership==1,1,0)
data$religion_peruana <- ifelse(data$church_membership==2,1,0)
data$religion_catholic <- ifelse(data$church_membership==3,1,0)
data$religion_maranatha <- ifelse(data$church_membership==4,1,0)
data$religion_israelita <- ifelse(data$church_membership==6,1,0)
data$religion_other <- ifelse(data$church_membership==5 | data$church_membership>6,1,0)

restricted_d <- data[which(data$treatment_catholic==1),]
column1 <- round(rbind(table(restricted_d$religion_catholic)[2], 100*table(restricted_d$religion_catholic)[2]/ length(restricted_d[,1]), table(restricted_d$religion_adventist)[2], 100*table(restricted_d$religion_adventist)[2]/ length(restricted_d[,1]), table(restricted_d$religion_peruana)[2], 100*table(restricted_d$religion_peruana)[2]/ length(restricted_d[,1]), table(restricted_d$religion_maranatha)[2], 100*table(restricted_d$religion_maranatha)[2]/ length(restricted_d[,1]), table(restricted_d$religion_israelita)[2], 100*table(restricted_d$religion_israelita)[2]/ length(restricted_d[,1]), table(restricted_d$religion_other)[2], 100*table(restricted_d$religion_other)[2]/ length(restricted_d[,1])),2)

restricted_d <- data[which(data$treatment_adventist==1),]
column2 <- round(rbind(table(restricted_d$religion_catholic)[2], 100*table(restricted_d$religion_catholic)[2]/ length(restricted_d[,1]), table(restricted_d$religion_adventist)[2], 100*table(restricted_d$religion_adventist)[2]/ length(restricted_d[,1]), table(restricted_d$religion_peruana)[2], 100*table(restricted_d$religion_peruana)[2]/ length(restricted_d[,1]), table(restricted_d$religion_maranatha)[2], 100*table(restricted_d$religion_maranatha)[2]/ length(restricted_d[,1]), table(restricted_d$religion_israelita)[2], 100*table(restricted_d$religion_israelita)[2]/ length(restricted_d[,1]), table(restricted_d$religion_other)[2], 100*table(restricted_d$religion_other)[2]/ length(restricted_d[,1])),2)

restricted_d <- data[which(data$treatment_mixed==1),]
column3 <- round(rbind(table(restricted_d$religion_catholic)[2], 100*table(restricted_d$religion_catholic)[2]/ length(restricted_d[,1]), table(restricted_d$religion_adventist)[2], 100*table(restricted_d$religion_adventist)[2]/ length(restricted_d[,1]), table(restricted_d$religion_peruana)[2], 100*table(restricted_d$religion_peruana)[2]/ length(restricted_d[,1]), table(restricted_d$religion_maranatha)[2], 100*table(restricted_d$religion_maranatha)[2]/ length(restricted_d[,1]), table(restricted_d$religion_israelita)[2], 100*table(restricted_d$religion_israelita)[2]/ length(restricted_d[,1]), table(restricted_d$religion_other)[2], 100*table(restricted_d$religion_other)[2]/ length(restricted_d[,1])),2)

restricted_d <- data[which(data$treatment_peruana==1),]
column4 <- round(rbind(table(restricted_d$religion_catholic)[2], 100*table(restricted_d$religion_catholic)[2]/ length(restricted_d[,1]), table(restricted_d$religion_adventist)[2], 100*table(restricted_d$religion_adventist)[2]/ length(restricted_d[,1]), table(restricted_d$religion_peruana)[2], 100*table(restricted_d$religion_peruana)[2]/ length(restricted_d[,1]), table(restricted_d$religion_maranatha)[2], 100*table(restricted_d$religion_maranatha)[2]/ length(restricted_d[,1]), table(restricted_d$religion_israelita)[2], 100*table(restricted_d$religion_israelita)[2]/ length(restricted_d[,1]), table(restricted_d$religion_other)[2], 100*table(restricted_d$religion_other)[2]/ length(restricted_d[,1])),2)

restricted_d <- data[which(data$treatment_maranatha==1),]
column5 <- round(rbind(table(restricted_d$religion_catholic)[2], 100*table(restricted_d$religion_catholic)[2]/ length(restricted_d[,1]), table(restricted_d$religion_adventist)[2], 100*table(restricted_d$religion_adventist)[2]/ length(restricted_d[,1]), table(restricted_d$religion_peruana)[2], 100*table(restricted_d$religion_peruana)[2]/ length(restricted_d[,1]), table(restricted_d$religion_maranatha)[2], 100*table(restricted_d$religion_maranatha)[2]/ length(restricted_d[,1]), table(restricted_d$religion_israelita)[2], 100*table(restricted_d$religion_israelita)[2]/ length(restricted_d[,1]), table(restricted_d$religion_other)[2], 100*table(restricted_d$religion_other)[2]/ length(restricted_d[,1])),2)


table <- cbind(column1, column2, column3, column4, column5)
column6 <- round(rbind(sum(table[1,], na.rm=T),100*sum(table[1,], na.rm=T)/length(data[,1]), sum(table[3,], na.rm=T),100*sum(table[3,], na.rm=T)/length(data[,1]), sum(table[5,], na.rm=T),100*sum(table[5,], na.rm=T)/length(data[,1]), sum(table[7,], na.rm=T),100*sum(table[7,], na.rm=T)/length(data[,1]), sum(table[9,], na.rm=T),100*sum(table[9,], na.rm=T)/length(data[,1]), sum(table[11,], na.rm=T),100*sum(table[11,], na.rm=T)/length(data[,1])),1)
table <- cbind(table, column6)
row13 <- t(rbind(table(data$treatment_catholic)[2], table(data$treatment_adventist)[2], table(data$treatment_mixed)[2], table(data$treatment_peruana)[2], table(data$treatment_maranatha)[2], length(data[,1])))
table <- round(rbind(table, row13),1)
colnames(table) <- c("Control", "Adventist", "Mixed", "Peruna", "Maranatha", "Total")
rownames(table) <- c("Catholic", "%", "Adventist", "%", "Peruana", "%", "Maranatha", "%", "Isralita", "%", "Other", "%", "Total")
table



